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" ^j" ■ We report on three topics on finite density simulations: (i) the derivative method for 

^— ^ hadronic quantities, (ii) phase fluctuations in the vicinity of the critical temperature and 

f—*) , (iii) the density of states method at finite isospin density. 

<n : 

^"V 81. Introduction 

Lattice QCD simulations at finite chemical potential li are extremely difficult due 

\ to the sign problem. Recently it has been realized that at small chemical potential 

one can study density effects on physical quantities by various approaches. x ) One 

. of the approaches is the derivative method which has been used for the study of 

response of meson masses with respect to li by QCD-TARO collaboration. 2 )' 3 ) The 
original idea of the derivative method may date back to the study of the fermion 
number susceptibilities 4 ) where the derivative of the fermion number density with 

\ repect to li were calculated. In section 2, we report on results from the derivative 

method for the response of meson masses as well as the chiral condensate. 

, At large li we believe that the standard Monte Cairo method based on impor- 

tance sampling fails due to the sign problem, i.e. the phase fluctuation is expected to 
be large for large li. However, we do not know exactly how the phase fluctuates with 
various simulation parameters. In section 3, we give results of the phase fluctuation 
in the vicinity of the critical temperature. We also present results of {tpip) an d the 
^3 ■ Polyakov loop for li < 0.25. 

J> , The most lattice simulations have been performed by methods based on impor- 

| tance sampling. There exists an alternative method called density of states method. 

An advantage of the density of states method is that one can obtain results for 



various values of simulation parameters without making independent simulations al- 
though the lattice size may be limited to a small one. In section 4, we give results 
from the density of states method for isospin density. Since the system at isospin 
density has no sign problem the results are compared with those from the standard 
Monte Carlo method. 



§2. Derivative method for hadronic quantities 



The derivative method extracts physical information at a small chemical poten- 
tial without directly simulating a system at finite chemical potential. The original 
idea of the derivative method comes from the calculations of the fermion number 
susceptibilities. 4 ) The QCD-TARO collaboration 2 )' 3 ) studied the response of meson 
masses by the derivative method. The method was also used for the study of the 
pressure. 5 )' 6 ) 
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Let us consider an observable O(ll). The expectation value of the observable 
O(fi) is given by 

fOAe'^dU , . 

<°>=-T-a ^-J7T> 2' 1 

/ Ae~ s sdU y ' 

where S g is the gluonic action and A stands for the fermion determinant. For two 
flavors of staggered quarks (u and d quarks), A is written as 

A = det M(pL u ) l l A det M{Li d ) 1/4 (2-2) 
= exp{^TrlnM(n u ) + jTrlnM(Li d )). (2-3) 

The first derivative of (O) with respect to \i at li = is given by 

where we used O = for simplicity. Note that here we used (^) = at li = 0. 
Similarly the second derivative at li = is given by 

d 2 < O > ■• • A A 

g^2 =<0 + 20->-<Oo- > cc , (2-5) 

where < Ao B > cc =< AB > - < A >< B >. 

Next let us consider the spatial hadronic correlator, 

C(x) = Y,(H(x,y,z,t)H(0,0,0,0)l}, (2-6) 

y,z,t 

and take derivatives of C(x). We assume that C(x) is dominated by a single pole 
contribution, 

C{x) = A(e~ Mx + e~ M{L - x) ), (2-7) 

where M is the hadron mass and L is the lattice size in the x direction. Taking the 
first and the second derivatives of C (x) with respect to li we obtain 

1 dC(x) 1 dA dM L. , ri ,. L L. , , 

" TTTTT + ^77T x {(* " T7)tanh[M(x - -)] - -}, (2-8) 



C(x) dii AdLi d/x 2' ' L v 2 n 2 

and 



1 d 2 C(x) Id 2 A f 2dAdM d 2 M x(l L , L.. L a . 

+ (tTTT^TT + "^tHO* " o ) tanh[M(x - -)] - -J2-9) 



C(x) dfx 2 AdLi 2 A djj, dji dii 2 2' 2 

+ (^) 2 {(- " \f + X " " |) tanh t M ^ " | 
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The left hand sides of these equations are used as fitting functions to the Monte 

Carlo data given by eq.(2-4) and (2-5). The response of masses with respect to jjl, 

dM , d 2 M , „ 

i.e. — — and are given as the fatting parameters. 

a\x dfi 2 

The QCD-TARO collaboration studied the response of the pseudo scalar (PS) 
meson screening mass for rif = 2 flavors on f 6 x 8 2 x 4 lattices. The first response of 
the PS meson mass turned out to be consistent with zero. Thus the first relevant re- 
sponse is the second one. Figure f shows the second responses of the PS meson mass. 
Here the isoscalar and isovector 
chemical potentials are defined by 

Vs = Vu = fid, (2-10) 

and 

Hv = Vu = -Vd, (2T1) 



-6 
-8 

















3E--- 


- - - at 












/ 






/ 






/ 

/ 






/ 










3E-' 


--at ^ |l v " 



Fig. 1. Second responses d 2 M/dfi 2 s and d 2 M/d(iy of 
the PS meson mass at m„a — 0.025 as a function 



respectively. The isovector chem- "a ~ 2 
ical potential is also called isospin -4 
chemical potential and we also use 
fly to stand for the isospin chem- 
ical potential. 

In the low temperature phase, - 10 080 o.ss 0.90 0.95 1.00 1.05 1.10 1.15 

the dependence of the mass on jig T/T c 
is small. This behavior is to be ex- 
pected since below T c the PS me- 
son is a Goldstone boson and per- of T/T c 2),3) 
sists its zero mass feature. On the 
other hand, above T c the PS me- 
son loses the Goldstone nature and can obtain its mass, as a result d 2 M/d[i 2 s seems 
to remain finite. 

The PS meson on [iy seems to behave differently. In the low temperature phase 
the PS meson mass tends to decrease with fiy. In the high temperature phase the 
PS meson also seems to decrease but the rate of the decrease is small. 

Next we consider the derivatives of {fiif)}. At \x = the first derivatives of (iptf)) 
with respect to both isoscalar and isovector chemical potentials are identically zero, 
i.e. 



dfi s dfi v 



= 0. (2-12) 



Thus the first relevant term is the second derivative. Figure 2 shows the second 

derivatives of {'tpip) with respect to us and nv f° r n f = 2 flavors on 16 x 8 2 x 4 

lattices, calculated by Choe et al. 2 ^' 7 ^ No significant difference can be seen between 

d 2 M) , d 2 {4>^) m , c , . , 

— —^5 — and — — 5 — . The strength of the responses increases as the quark mass 

decreases. The similar results are also obtained for the NJL model. 8 ) 
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Fig. 2. (left) and ^fi- (right) as a function of T/T c . 7) 



dp% 



dpi 



The results of the second derivatives can be used for evaluating (iptp) at small 
chemical potentials as 



l d 2 W) 2 
2 V M ' 



(2-13) 



Figure 3 shows (ipip(ns)) extrapolated to finite fis using eq.(2-13). We see that 
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Fig. 3. (ipip) extrapolated to finite ps at m q = 0.025. 7 - 1 Here T c is the critical temperature at 
p = 0. 



qualitatively the critical temperature decreases as fi increases. Similar behavior can 
be seen for the case of the isovector chemical potential. 



5 



§3. Phase fluctuation in the vicinity of the critical temperature 

There is no satisfactory method to simulate the system at large /i. In principle 
the reweighting method can be applied for any /x as 



<o(m)> 



<0(M)e* 



(e ie ) 



(3-1) 



AM' 



where (0)a v stands for the expectation value of the operator O in an ensemble 
at isospin density, i.e. the configurations are generated with the phase quenching 
measure ~ \det M(^i)\ n f/ 4 e~ s 9 . For large /i eq.(3-l) is expected to be impractical 
since {e ie )^ v becomes small. For such small (e l6 )/j, v , in order to obtain a meaningful 
value of eq.(3T) one needs extremely high statistics. Furthermore there is another 
difficulty: the calculation of the complex phase 9 contains a determinant calculation 
which is computationally costly. For small \i one may use the Taylor expansion of 
the determinant as 9 ) 



detM(/i) = exp(TWnM(0) + Tr 



1 dM 



M d[i 



(3-2) 



Ai=0 



The first term of eq.(3-2) is real and the second term is pure imaginary. Thus 9 at 
small fi is given by 

9 = 6 + 0(^), (3-3) 



where 



6 = ImTr 



1 dM 



M d\x 



x ii. 



(3-4) 



In this expression no determinant calculation is required. Instead, 9 is given by a 
trace calculation and the computational cost is much reduced. An empirical study 
shows that the quality of the approximation 9 is valid for /x < O.2. 10 ) Figure 4 shows 9 



phase 



vs the exact 

The agreement between 9 and 9 
is excellent for fj, = 0.1 and 0.2. 3 
However for \x = 0.3, 9 deviates 
from the exact value 9. 2 
For large /i, in order to ob- 1 
tain 9, using the approximation J 
is not sufficient and one needs -3 
the determinant calculation to ob- * _ 1 
tain the value of 9. Sasai, Naka- 
mura and Takaishi studied the " 2 
phase fluctuations by calculating 9 _ 3 
without any approximation. Fig- 
ure 5 shows the phase fluctuation 
(cos(# /2)) fj, v for nj = 2 as a func- Fig. 4. 6 vs the exact phase 9 for 3 values of fj, 

tion of fi. n ) (cos{9/2))^ v is the 



9. 




mu=0.1 

mu=0.2 
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10) 
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average value over the phase quenched ensemble. The phase fluctuation increases, 
i.e., (cos(0/2)) MV decreases as /i increases. However in the deconfinement region, the 
phase fluctuation is much smaller than that in the confinement region. 

SU(3) N=2 m=0.05 8 3 x4 lattice 

I i i ' V 1 S3 

i 

i i 
i 

M 



x n=0.1 

+ ^=0.2 : 

p.=0.25 ; 

i i i I i i i i I i i i i I i i i i L 

"5.2 5.25 5.3 5.35 5.4 
P 

Fig. 5. Phase fluctuation (cos(0/2)) MV for n/ = 2 in the isospin ensemble. 11 ' /3 c (m = 0) ~ 5.32. 

With the simulation parameters used here, the phase fluctuation is not signifi- 
cantly large and the expectation values at finite /x can be obtained by the reweighting 
from the isospin ensemble as eq.(3-l). Figure 6 shows (ipip) and Polyakov loop with 
and without reweighting for /i < 0.25. We do not see any difference between the 
results with and without reweighting. This is consistent with the result of {$ip) in 
12). This might indicate that for [i < 0.25 the phase effects are small. This fact is 
also consistent with that the phase diagram obtained at small [xy is similar to that at 
small /is. 13 )' 14 ) The phase diagram at small \iy was calculated for n/ = 2 staggered 
fermions 14 ) and also for the Wilson fermions 13 ) with the plaquette and DBW2 15 )~ 17 ) 
gauge actions. 

Although the isospin system has a positive measure and can be simulated with 
the standard Monte Carlo methods, for large chemical potentials we face the com- 
putational difficulty: the matrix solver does not converge. As a result the simulation 
can not be performed. This feature is also mentioned in 12). Interestingly for fur- 
ther large chemical potentials /u > /x c (/u c is a certain value which is dependent of 
the simulation parameters), the matrix solver converges. 18 ) 

§4. Density of states method at isospin density 

The most lattice simulations are performed by the Monte Carlo methods based 
on importance sampling. Here we report on an alternative approach, the density of 
states (DOS) method. We apply the DOS method for isospin density and compare 
the results with those from the standard Monte Carlo method. 
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Fig. 6. (ipip) (left) and Polyakov loop (right) with and without reweighting. 11 - ) "no phase" stands 
for the results without reweighting. 



The DOS method has been applied for gauge theories 19 ^ and QED with dy- 
namical fermions. 20 ) Luo 21 ^ applied the DOS method for QCD and argued that if 
the eigenvalues of the Dirac operators are determined the thermodynamic quantities 
derived from the partition function can be evaluated at any quark mass and flavor. 

The expectation value of the operator O is given by 

(0) = | J[dU)0[U]detM(m q ,^ n f^eM-/3S g [U)), (4-1) 

where M(m q , fi) is assumed to be a staggered fermion matrix at quark mass m q and 
at chemical potential \x. We define n(E) as 



n (E) = J [dU}8(6VE - S g [U\), (4-2) 

where V is the number of lattice sites and E is the plaquette energy. 
Using n(E), (O) is expressed by 

(O) = ^- f dEn(E)e- l36VE (O det M{fi) N f/ 4 ) E , (4-3) 
Z n J 

Z n = J dEn{E)e- (36VE (detM^) N f/ 4 ) E . (4-4) 

(•)e stands for the microcanonical averages with fixed E. If these microcanonical 
averages are determined as a function of E, (O) is given at any f3. Furthermore, Luo 
argued that if one stores the eigenvalues of for all configurations, then one can 



where 
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evaluate (det M(fj,) N f^)E at any quark mass and flavor as 



N C V 



(detM(ti) Nf/i }E 



(4-5) 



where Aj(//) is the i-th eigenvalue of the massless fermion matrix M{m q = 0) and 
iV c = 3 for SU(3). On the other hand, in general (O det M(u) N fl A ) jg is not calculable 
at any quark mass and flavor. However for O = ipip(fi) which is given with the trace 
of M _1 one can also evaluate it at any quark mass and flavor as 



((^TrMi^detM^^E = E 



N C V 



Xi(li) + m q 



(II (W + 



q )) Nf/A )E. 

(4-6) 



Using these equations one can obtain (i[)ij){u)) at any (3, quark mass and flavor. 

Here we comment on the con- 
struction of the density of states. 
One can define the DOS for any 
other quantities. For instance, 
Gocksch 22 ) constructed the DOS for 
the complex phase. Ambjorn et 
a/. 23 - 1 constructed the DOS for the 
number density in their factorization 
method. In these definitions of the 
DOS, simulation parameters /?, ray 
and m q are absorbed in the DOS 
and the simulation parameters are 
not variable. 

In the following we show results 
on 4 4 lattices at isospin chemical 
potential [iy . The DOS n(E) in 
eq.(4-2) can be obtained using the 
quenched data as 21 ) 




Fig. 7. 
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(4-7) 



Figure 7 shows — In n(E) /V as a function of plaquette energy E. The time consuming 
part of the method is the calculations of (det M(n) N f/ 4 ) E and (O det M(n) N f/ 4 ) E 
which contain the eigenvalue calculation. In order to generate configurations at fixed 
E we used the over-relaxation method. At each E we generated 100 configurations. 
Each configuration is separated by 100 over-relaxed updates. For each configuration 
we calculate eigenvalues and those eigenvalues are used to evaluate eqs.(4-5) and 
(4-6). Figure 8 shows (det M(ny) N f^)E at uy = 0.2 and m q = 0.025 for various 
flavors as examples. 

The isospin system has a positive measure and can be simulated with the stan- 
dard Monte Carlo algorithm as R-algorithm. 25 ) Figure 9 compares the results from 
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the DOS method with those from the R-algorithm. They are in good agreement 
with each other. Only a small difference is seen in the phase transition region. 
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Fig. 9. (ipip) for rif = 2 at \iv = 0.2 and m q — 0.05 



Figure 10 shows (ijjip) for different nj at /xy = 0.25 and m q = 0.05. One can see 
that how p c changes as n/, i.e. f3 c decreases as nj increases. 

In the DOS method one may take various combinations of parameters. Let us 
consider the case of nj = 1 + 1 with non-degenerate quark masses m\ and mi- In 
this case one must calculate the following microcanonical averages: 

(| det A( mi )\ N f/ 4 \ det A(m 2 )\ N f/ 4 ) E , (4-8) 

(#(m i=li2 )| det Z\(mi)|^/ 4 | det A(m 2 ) \ N i' A ) E . (4-9) 

Since the eigenvalues are stored, it is easy to calculate these microcanonical averages. 
On the other hand, in the conventional algorithm as R-algorithm, it is not easy to 
simulate the non-degenerate system since one needs a different program to perform 
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Fig. 10. {ip4>} f° r different rif at fiv = 0.25 and m q = 0.05 
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the simulation. Figure ll(left) shows (if^ip) for n/ = 1 + 1 with different quark masses 
( mi = 0.05 and m 2 = 0.025 ) at fi v = 0.2. 

Similarly one can also consider non-degenerate isospin chemical potentials m 
and fi2- In this case one calculates 

(Idet^^Or/^ldet^^)!^/ 4 )^, (4-10) 

(^(^=i, 2 )|detZ\( / x 1 )| Af // 4 |det^(/i 2 )r //4 )i?. (4-11) 

Figure 11 (right) shows (tpip) for nj = 1 + 1 with different isospin chemical potentials 
( fix = 0.2 and /x 2 = 0.3 ) at m q = 0.025. 
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Fig. 11. {if>i{j(mq, [iv)) for Nf = 1 + 1 at \iv = 0.2 with different quark masses, (m q = 0.05 and 



0.025) (left), and at m q = 0.025 with different chemical potentials (^=0.2 and 0.3) 
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In the DOS method one can easily obtain results for various parameters without 
making independent simulations. Thus the DOS method is considered to be useful 
to explore a wide parameter space. On the other hand, the eigenvalue calculations 
are computationally difficult on large lattices. Therefore the application of the DOS 
method might be limited on small lattices. 
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